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We investigate dissipative phase transitions in an open central spin system. In our model the 
central spin interacts coherently with the surrounding many-particle spin environment and is subject 
to coherent driving and dissipation. We develop analytical tools based on a self-consistent Holstein- 
Primakoff approximation that enable us to determine the complete phase diagram associated with 
the steady states of this system. It includes first and second-order phase transitions, as well as 
regions of bistability, spin squeezing and altered spin pumping dynamics. Prospects of observing 
these phenomena in systems such as electron spins in quantum dots or NV centers coupled to lattice 
nuclear spins are briefly discussed. 



I. INTRODUCTION 

Statistical Mechanics classifies phases of a given system 
in thermal equilibrium according to its physical prop- 
erties. It also explains how changes in the system pa- 
rameters allow us to transform one phase into another, 
sometimes abruptly, which results in the phenomenon of 
phase transitions. A special kind of phase transitions oc- 
curs at zero temperature: such transitions are driven by 
quantum fluctuations instead of thermal ones and are re- 
sponsible for the appearance of exotic quantum phases in 
many areas of physics. These quantum phase transitions 
have been a subject of intense research in the last thirty 
years, and are expected not only to explain interesting 
behavior of systems at low temperature, but also to lead 
to new states of matter with desired properties (e.g., su- 
perconductors, -fluids and -solids, topological insulators 

mm)- 

Phase transitions can also occur in systems away from 
their thermal equilibrium. For example, this is the case 
when the system interacts with an environment and, 
at the same time, is driven by some external coherent 
source. Due to dissipation, the environment drives the 
system to a steady state, po(g), which depends on the 
system and environment parameters, g. As g is changed, 
a sudden change in the system properties may occur, giv- 
ing rise to a so called dissipative phase transition (DPT) 
[7HT2] . DPT have been much less studied than traditional 
or quantum ones. With the advent of new experimental 
techniques that allow to observe them experimentally, 
they are starting to play an important role [13] . Moreover 
they offer the intriguing possibility of observing critical 
effects non-destructively because of the constant intrin- 
sic exchange between system and environment |14j . In 
equilibrium statistical mechanics a large variety of toy 
models exist that describe different kind of transitions. 
Their study lead to a deep understanding of many of 
them. In contrast, in the case of DPT few models have 
been developed. 

The textbook example of a DPT occurs in the Dicke 



model of resonance fluorescence [3 [15] . There, a system 
of spins interacts with a thermal reservoir and is exter- 
nally driven. Experimental [IB] and theoretical studies 
[P7H20] revealed interesting features such as optical mul- 
tistability, first and second order phase transitions, and 
bipartite entanglement. 

In this paper, we analyze another prototypical open 
system: The model is closely related to the central spin 
system (CSS) which has been thoroughly studied in ther- 
mal equilibrium [21-23 . In its simplest form, it consists 
of a set of spin-1/2 particles (in the following referred to 
as the nuclear spins), uniformly coupled to a single spin- 
1/2 (referred to as the electron spin). In the model we 
consider, the central spin is externally driven and decays 
through interaction with a Markovian environment. Re- 
cently, the CSS model has found application in the study 
of solid state systems such as electron and nuclear spins 
in a quantum dot [23j or an Nitrogen- Vacancy-center. 

In what follows we first provide a general framework 
for analyzing DPT in open systems. In analogy with the 
analysis of low energy excitations for closed systems, it is 
based on the study of the excitation gap of the system's 
Liouville operator C. We illustrate these considerations 
using the central spin model. For a fixed dissipation 
strength 7, there are two external parameters one can 
vary, the Rabi frequency of the external driving field, 
51, and the Zeemann shift, u. We present a complete 
phase diagram as a function of those parameters, char- 
acterize all the phases, and analyze the phase transitions 
occurring among them. To this end, we develop a series 
of analytical tools, based on a self-consistent Holstein- 
Primakoff approximation, which allows us to understand 
most of the phase diagram. In addition, we use numerical 
methods to investigate regions of the diagram where the 
theory yields incomplete results. Combining these tech- 
niques, we can identify two different types of phase tran- 
sitions, and regions of bistability, spin squeezing, and en- 
hanced spin polarization dynamics. We will also identify 
regions where anomalous behavior occurs in the approach 
to the steady state. Intriguingly, recent experiments with 
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quantum dots, in which the central (electronic) spin is 
driven by a laser and undergoes spontaneous decay, re- 
alize a situation very close to the one we study here and 
show effects such as bistability, enhanced fluctuations, 
and abrupt changes in polarization in dependence of the 
system parameters [211 [25] . 

This paper is organized as follows. Section [II] sets the 
general theoretical framework underlying our study of 
DPT. Section |III| introduces the model, and contains a 
structured summary of the main results. In Section IV 
we develop the theoretical techniques and use those tech- 
niques to analyze the various phases and classify the dif- 
ferent transitions. Thereafter in Section [V] numerical 
techniques are employed to explain the features of the 
phase diagram which are not captured by the previous 
theory. Possible experimental realizations and a general- 
ization of the model to inhomogeneous coupling are dis- 



cussed in Section VI Finally we summarize the results 
and discuss potential applications in Section |VII[ 



II. GENERAL THEORETICAL FRAMEWORK 

The theory of quantum phase transitions in closed sys- 
tems is a well established and extensively studied area in 
the field of statistical mechanics. The typical scenario 
is the following: a system is described by a Hamilto- 
nian, H(g), where g denotes a set of systems parame- 
ters (like magnetic fields, interactions strengths, etc). At 
zero temperature and for a fixed set of parameters, g, the 
system is described by a quantum state, i[>o(g), fulfilling 
[H(g) - E^ Q (g)}\i> Q {g)) = 0, where E^{g) is the ground 
state energy. As long as the Hamiltonian is gapped (i.e., 
the difference between Eo(g) and the first excitation en- 
ergy is finite), any small change in g will alter the physical 
properties related to the state \tpo(g)) smoothly and we 
remain in the same phase. However, if the first excitation 
gap A = E l f, 1 (g) — E^ a (g) closes at a given value of the 
parameters, g = go, it may happen that the properties 
change abruptly, in which case a phase transition occurs. 

In the following we adapt analogous notions to the case 
of DPT and introduce the concepts required for the sub- 
sequent study of a particular example of a generic DPT 
in a central spin model. 

We consider a Markovian open system, whose evolu- 
tion is governed by a time-independent master equation 
p = C(g)p. The dynamics describing the system are con- 
tractive implying the existence of a steady state. This 
steady state po(g) is a zero eigenvector to the Liouville 
superoperator C(g)po(g) = 0. This way of thinking par- 
allels that of quantum phase transitions, if one replaces 
[H(g) — E^ (g)] — > C(g). Despite the fact that these 
mathematical objects are very different (the first is a Her- 
mitian operator, and the second a hermiticity-preserving 
superoperator) , one can draw certain similarities between 
them. For instance, for an abrupt change of pa(g) (and 
thus of certain system observables) it is necessary that 
the gap in the (in general complex) excitation spectrum 



of the system's Liouville operator C(g) closes. The rele- 
vant gap in this context is determined by the eigenvalue 
with largest real part different from zero (it can be shown 
that Re(A) < for all eigenvalues of £ [IS])- The van- 
ishing of the real part of this eigenvalue - from here on 
referred to as asymptotic decay rate (ADR) [27] - in- 
dicates the possibility of a non-analytical change in the 
steady state and thus is a necessary condition for a phase 
transition to occur. 

In our model system, the Liouvillian low excitation 
spectrum, and the ADR in particular, can in large parts 
of the phase diagram be understood from the complex 
energies of a stable Gaussian mode of the nuclear field. 
We find first order transitions where the eigenvalue of 
this stable mode crosses the eigenvalue of a metastable 
mode at zero in the projection onto the real axis. The 
real part of the Liouvillian spectrum closes directly as the 
stable mode turns metastable and vice versa. A finite dif- 
ference in the imaginary parts of the eigenvalues across 
the transition prevents a mixing of the two modes and 
the emergence of critical phenomena such as a change in 
the nature of the correlations in the steady state at the 
critical point. In contrast, we also find a second order 
phase transition where the ADR vanishes asymptotically 
as both mode energies become degenerate (at zero) in the 
thermodynamic limit. Mixing of the two modes at the 
critical point gives rise to diverging correlations in the nu- 
clear system. This observation parallels the classification 
of quantum phase transitions in closed systems. There, 
a direct crossing of the ground and first excited state en- 
ergy for finite systems (mostly arising from a symmetry 
in the system) typically gives rise to a first order phase 
transition. An asymptotical closing of the first excita- 
tion gap of the Hamiltonian in the thermodynamic limit 
represents the generic case of a second order transition 

Besides the analogies described so far [cf. Table [I] , 
there are obvious differences, like the fact that in DTP 
Po (g) may be pure or mixed, and that some of the char- 
acteristic behavior of a phase may also be reflected in 
how the steady state is approached. Non-analyticities 
in the higher excitation spectrum of the Liouvillian are 
associated to such dynamical phases. 



III. MODEL AND PHASE DIAGRAM 

A. The Model 

We investigate the steady state properties of a homoge- 
neous central spin model. The central spin - also referred 
to as electronic spin in the following - is driven resonantly 
via suitable optical or magnetic fields. Dissipation causes 
electronic spin transitions from the spin-up to the spin- 
down state. It can be introduced via standard optical 
pumping techniques. [291 130]. Furthermore, the central 
spin is assumed to interact with an ensemble of ancilla 
spins - also referred to as nuclear spins in view of the 
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[H-Ei D ]\ih) = 


po = argmin[||£p|| tr ] 
l|p|| tr =i 

Cp = 


Phase transition 


Non-analyticity in F(px) 


A = E^ 11 — E^, vanishes 


ADR = max[Re(A p )] vanishes 



Table I: Non-exhaustive comparison of thermal phase transitions (TPT), QPT and DPT. The concepts for DPT parallel in 
many respects the considerations for QPT and TPT. | • || tr denotes the trace norm and S the entropy. Note that if the steady 
state is not unique, additional steady states may come with a non-zero imaginary part of the eigenvalue and then appear in 
pairs: Cp = ±iyp (y £ K). 



mentioned implementations |23j- by an isotropic and ho- 
mogeneous Heisenberg interaction. In general this hyper- 
fine interaction is assumed to be detuned. Weak nuclear 
magnetic dipole-dipole interactions are neglected. 

After a suitable transformation which renders the 
Hamiltonian time-independent, the system under consid- 
eration is governed by the master equation 

p = C P (1) 
- Jj(S-pS + - ^{S+S-, P }) - i[H s + Hj + Hsi, p], 

where {•, •} denotes the anticommutator and 

H s = ,m(S+ + S-), (2) 
H 1 = Sul s , (3) 
H S i = a/2(S + r + S-I+) + aS + S-I z . (4) 

S a and I a = of ( a = +i — : z ) denote electron 
and collective nuclear spin operators, respectively. Jfl 
is the Rabi frequency of the resonant external driving 
of the electron (in rotating wave approximation), while 
5lu = lu — a/2 is the difference of hyperfine detuning lu 
and half the individual hyperfine coupling strength a. 8lu 
for instance can be tuned via a static magnetic fields in 
z direction. Note that Hj + Hsi = aSI + ujI z , describing 
the isotropic hyperfine interaction and its detuning. The 
rescaling of the electron driving and dissipation in terms 
of the total (nuclear) spin quantum number J |51j is in- 
troduced here for convenience and will be justified later. 
Potential detunings of the electron driving - correspond- 
ing to a term AS Z in the Hamiltonian part of the master 
equation - can be neglected if A <C Ja. 

In the limit of strong dissipation ]»a the electron 
degrees of freedom can be eliminated and Eq. ([!]) reduces 
to 

& := Tr s (p) = 7eff (/-a7+ - l -{I+r,a}) (5) 
-i [Q e ffly + SloI z ] , 



where 7 c ff = — , f2 e flf = and o is the reduced density 
matrix of the nuclear system. This is a generalization of 
the Dicke model of resonance fluorescence as discussed in 

nana can. 

Master Eq. ([I]) has been theoretically shown to display 
cooperative nuclear effects such as superradiance even for 
inhomogeneous electron nuclear coupling [31 . In analogy 
to the field of cooperative resonance fluorescence, the sys- 
tem's rich steady state behavior comprises various critical 
effects such as first and second order DPT and bistabil- 
ities. In the following we provide a qualitative summary 
of the phase diagram and of the techniques developed to 
study the various phases and transitions. 



B. Phenomenological Description of the Phase 
Diagram 

For a fixed dissipation rate 7 = [S5] the differ- 
ent phases and transitions of the system are displayed 
schematically in Fig. [l] in dependence on the external 
driving Q and the hyperfine detuning w. We concentrate 
our studies on the quadrant Q, u > in which all in- 
teresting features can be observed. In the following, we 
outline the key features of the phase diagram. 

First we consider the system along the line segment x 
(lo = u)q, < Clo), where Q,q = ojq = a/2 define a critical 
driving strength and critical hyperfine detuning, respec- 
tively. Here Hj vanishes and the steady state can be con- 
structed analytically as a zero-entropy factorized state of 
the electron and nuclear system. The nuclear field builds 
up to compensate for the external driving - forcing the 
electron in its dark state \l) - until the maximal polariza- 
tion is reached at the critical value 51o- Above this point 
the nuclear system cannot compensate for the driving f2 
anymore and a solution of different nature, featuring fi- 
nite electron inversion and entropy is found. The point 
S7o features diverging spin entanglement and will be iden- 
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Figure 1: Schematic of the different phases and transitions 
of Master Eq. dTp . fn the two main phases of the system A 
(yellow) and B (blue) - which together cover the whole phase 
diagram - the system is found in a RSTSS (cf. text). While 
phase A is characterized by normal spin pumping behavior 
(large nuclear polarization in the direction of the dissipation) 
and a low effective temperature, phase B displays anomalous 
spin pumping behavior (large nuclear polarization in opposing 
direction to the dissipation) and high temperature. They are 
separated by the first order phase boundary b which is asso- 
ciated with a region of bistability C (framed by the boundary 
c). Here a second non-Gaussian solution appears, besides the 
normal spin pumping mode of A. The region of bistability C 
culminates in a second order phase transition at (wo, fio)- Be- 
low this critical point the system is supercritical and no clear 
distinction between phases A and B exists. In this region a 
dynamical phase D emerges, characterized by anomalous be- 
havior in the approach to the steady state. For a detailed 
description of the different phases and transitions see Sec- 
tion utm 



tified as a second order phase transition. 

For the separable density matrix p n = (01, \ip) — 
\\) <X> | a) the only term in Master Eq. which is not 
trivially zero is the Hamiltonian term S + (^I~ + JO). 
However, choosing |a) as an approximate eigenstate of 
the lowering operator I~ \a) « a \a) (up to second order 
in e = l/VJ) with a = — 2JQ/a = — JJ7/f2o (o is the in- 
dividual hyperfine coupling constant) the corresponding 
term in E q. Qip vanishes in the thermodynamic limit. In 
Appendix |A 1| we demonstrate that approximate eigen- 
states \a) can be constructed as squeezed and displaced 
vacua in a Holstein-Primakoff |32j picture up to a cor- 
rection of order 1/3. The squeezing of the nuclear state 
depends uniquely on the displacement such that these 
states represent a subclass of squeezed coherent atomic 
states |33| . Remarkably, this solution - where along the 
whole segment x the system settles in a separable pure 
state - exists for all values of the dissipation strength 7. 



In the limit of vanishing driving 51 = the steady state 
trivially is given by the fully polarized state (being the 
zero eigenstate of the lowering operator), as the model 
realizes a standard optical spin pumping setting for dy- 
namical nuclear polarization |34j . With increasing f2, the 
collective nuclear spin is rotated around the y-axis on the 
surface of the Bloch sphere such that the effective Over- 
hauser field in x-direction compensates exactly for the 
external driving field on the electron spin. As a conse- 
quence along the whole segment x the dissipation forces 
the electron in its dark state ||), and all electron observ- 
ables, but also the entropy and some nuclear observables, 
are independent of SI. 

Furthermore, the steady state displays increased nu- 
clear spin squeezing in y-direction (orthogonal to the 
mean polarization vector) when approaching the critical 
point. A common measure of squeezing is defined via the 
spin fluctuations orthogonal to the mean polarization of 
the spin system. A state of a spin- J system is called spin 
squeezed |33j . if there exists a direction n, orthogonal to 
the mean spin polarization (/), such that: 



£i = 2<AJ3)/|$|<l. 



(6) 



In [35] it was shown that every squeezed state also con- 
tains entanglement among the individual constituents. 
Moreover, if ^ < ^ then the spin squeezed state con 



tains A:-particle entanglement [351 [37J . In Appendix A 1 
we show that the squeezing parameter in j/-direction 
for an approximate I~ eigenstate \a) is given as £| = 

y/1 - a 2 / J 2 +0(1/ J) = y/l - (n/n ) 2 +O(l/J). Note 
however, that this equation is valid only for £| > \j\fj . 
For higher squeezing the operator expectation values con- 
stituting the term of order 0(1/ J) can attain macro- 
scopic values of order \/J . For SI < S7o we find that 
the nuclear spins are in a highly squeezed minimum un- 
certainty state, with fc-particle entanglement [53 . Close 
to the critical point k becomes of the order of \fj 
[£§ = 0(1 /\fj)} indicating diverging entanglement in 
the system. 

Since the lowering operator is bounded (||/ _ || < J), 
at Q = fio where the nuclear field has reached its maxi- 
mum value, the zero entropy solution constructed above 
ceases to exist. For large electron driving, where £1 3> ^0 
sets the dominant energy scale, the dissipation 7 results 
in an undirected diffusion in the dressed state picture 
and in the limit 17 — > 00 the system's steady state is 
fully mixed. In order to describe the system for driv- 
ing strength > fi , in Section [IV A| we develop a per- 
turbative theory designed to efficiently describe a class 
of steady states where the electron and nuclear spins 
are largely decoupled and the nuclear system is found 
in a fully polarized and rotated state with, potentially 
squeezed, thermal Gaussian fluctuations (also referred 
to as rotated squeezed thermal spin states - RSTSS or 
the Gaussian mode) It is fully characterized by its mean 
polarization as well as the spin squeezing and effective 
temperature T e g of the fluctuations (cf. Appendix [B| . 
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Squeezed coherent atomic states, which constitute the 
solution along segment x, appear as a limiting case of 
this class for zero temperature T ef j = 0. 

We conduct a systematic expansion of the system's Li- 
ouville operator in orders of the system size l/VJ, by 
approximating nuclear operators by their semiclassical 
values and incorporating bosonic fluctuations up to sec- 
ond order in an Holstcin-Primakoff picture. The result- 
ing separation of timescales between electron and nuclear 
dynamics is exploited in a formalized adiabatic elimina- 
tion of the electron degrees of freedom. The semiclassical 
displacements (i.e., the electron and nuclear direction of 
polarization) are found self-consistently by imposing first 
order stability of the nuclear fluctuations. For a given 
set of semiclassical solutions we derive a second order re- 
duced master equation for the nuclear fluctuations which, 
in the thermodynamic limit, contains all information on 
the nuclear state's stability, its steady state quantum 
fluctuations and entanglement as well as the low exci- 
tation dynamics in the vicinity of the steady state and 
thus allows for a detailed classification of the different 
phases and transitions. 

Using this formalism, we find that the system enters a 
new phase at the critical point Qq, in which the nuclear 
field can no longer compensate for the external driving, 
leading to a finite electron inversion and a nuclear state 
of rising temperature for increasing driving strength. At 
the transition between the two phases, the properties 
of the steady state change non-analytically and in Sec- 
tion IV B 2| we will find an asymptotic closing of the Li- 
ouvillian gap (cf. Section [TTJ) at the critical point, as the 
Liouvillian's spectrum becomes continuous in the ther- 
modynamic limit. We will characterize the critical point 
(cjo, Oo) as a second order phase transition. 

Allowing for arbitrary hypcrfinc detunings w, a phase 
boundary emerges from the second order critical point 
(line b in Fig. [T]) , separating two distinct phases A (blue) 
and B (red) of the Gaussian mode. The subregion C of 
A indicates a region of bistability associated to the phase 
boundary b and is discussed below. 

At 51 = the semiclassical equations of motion feature 
two steady state solutions. Not only the trivial steady 
state of the spin pumping dynamics - the fully polar- 
ized state in — z direction - but also an inverted state 
where the nuclear system is fully polarized in +z direc- 
tion is a (unstable) solution of the semiclassical system. 
Quantum fluctuations account for the decay of the lat- 
ter solution of anomalous spin pumping behavior. The 
two semiclassical solutions (the corresponding quantum 
states are from here on referred to as the normal and 
anomalous spin pumping mode, respectively) persist for 
finite f2. As we show employing the formalism described 
above (Section IV B 3), quantum fluctuations destabilize 
the mode of anomalous behavior in region A of the phase 
diagram. The stable Gaussian solution in phase A dis- 
plays a behavior characterized by the competition of dis- 
sipation 7 and the onsetting driving field Q. The nu- 
clear state is highly polarized in the direction set by the 



decay, and the electron spin starts aligning with the in- 
creasing external driving field. Furthermore, the normal 
spin pumping mode of phase A is characterized by a low 
effective spin temperature. 

The analysis of the low excitation spectrum of the Li- 



ouvillian (Section IV B 4 1 shows a direct vanishing of the 



ADR at the phase boundary b between A and B, while 
the imaginary part of the spectrum is gapped at all times. 
At this boundary, the normal mode of phase A destabi- 
lizes while at the same the metastable anomalous mode 
turns stable defining the second phase B. The two mode 
energies are non-degenerate across the transition prevent- 
ing a mixing of the two modes and the emergence of crit- 
ical phenomena such as diverging entanglement in the 
system. Phase B - anomalous spin pumping - is char- 
acterized by a large nuclear population inversion, as the 
nuclear field builds up in opposite direction of the dis- 
sipation. At the same time the electron spin counter 
aligns with the external driving field Q. In contrast to 
the normal mode of phase A, phase B features large fluc- 
tuations (i.e., high effective temperature) in the nuclear 
state, which increase for high f2, until at some point the 
perturbative description in terms of RSTSS breaks down 
and the system approaches the fully mixed state. Note 
that region A also transforms continuously to B via the 
lower two quadrants of the phase diagram (Fig. [T]) . In 
this supercritical region ^38 no clear distinction of the 
two phases exist. 

To complete the phase diagram, we employ numeri- 
cal techniques in order to study steady state solutions 
that go beyond a RSTSS description in Section |Vj The 
subregion of A labeled C indicates a region of bistability 
where a second steady state solution (besides the nor- 
mal spin pumping Gaussian solution described above) 
appears, featuring a non-Gaussian character with large 
fluctuations of order J. Since this mode cannot be de- 
scribed by the perturbative formalism developed in Sec- 
tion IV (which by construction is only suited for low fluc- 



tuations <C J) we use numerical methods to study this 
mode in Section [V] for finite systems. We find that the 
non-Gaussian mode (in contrast to the Gaussian mode 
of region A) is polarized in +z direction and features 
large fluctuations of the order of J. Additionally this 
solution displays large electron-nuclear connected corre- 
lations (Silj) — (Si)(Ij). It emerges from the anomalous 
spin pumping mode coming from region B and the sys- 
tem shows hysteretic behavior in region C closely related 
to the phenomenon of optical bistability [39] . 

A fourth region is found in the lower half of the phase 
diagram (D). In contrast to the previous regions, area 
D has no effects on steady state properties. Instead the 
region is characterized by an anomalous behavior in the 
low excitation dynamics of the system. The elementary 
excitations in region D are overdamped. Perturbing the 
system from its steady state, leads to a non-oscillating 
exponential return. This behavior is discussed at the 
end of Section |IV B 3| where we study the low excita- 
tion spectrum of the Liouvillian in this region within the 
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perturbative approach. 

In summary, all the phases and transitions of the sys- 
tem are displayed in Fig. [T] Across the whole phase di- 
agram one solution can be described as a RSTSS - a 
largely factorized electron-nuclear state with rotated nu- 
clear polarization and Gaussian fluctuations. Phase A 
hereby represents a region of normal spin pumping be- 
havior. The system is found in a cold Gaussian state, 
where the nuclear spins are highly polarized in the di- 
rection set by the electron dissipation and the electron 
spin aligns with the external driving for increasing field 
strength. In contrast, phase B displays anomalous spin 
pumping behavior. The nuclear system displays popula- 
tion inversion (i.e., a polarization opposing the electron 
pumping direction) while the electron aligns in opposite 
direction of the driving field. Furthermore the state be- 
comes increasingly noisy, quantified by a large effective 
temperature, which results in a fully mixed state in the 
limit of large driving strength £1 — > oo. Along segment x 
the state becomes pure and factorizes exactly with a nu- 
clear field that cancels the external driving exactly. The 
nuclear state can be described using approximate eigen- 
states of the lowering operator I~ which display diverg- 
ing squeezing approaching the second order critical point 
il . From this critical point a first order phase boundary 
emerges separating phases A and B. It is associated with 
a region of bistability (area C) , where a second solution 
appears featuring a highly non-Gaussian character. The 
system shows hysteretic behavior in this region. Region 
D is a phase characterized by its dynamical properties. 
The system shows an overdamping behavior approaching 
the steady state, which can be inferred from the excita- 
tion spectrum of the Liouvillian. 

Let us now describe the phases and transitions involv- 
ing the Gaussian mode in detail. 



include criticality in both steady state and low excita- 
tion spectrum, spin squeezing and entanglement as well 
as altered spin pumping dynamics. Whenever feasible 
we compare the perturbative results with exact diago- 
nalization techniques for finite systems and find excellent 
agreement even for systems of a few hundred spins only. 
First in Section pV B 2| we apply the developed theory ex- 
emplarily along the segment x, to obtain further insights 
in the associated transition at Slo- In Section |"IV B 3| we 
then give a detailed description of the different phases 
that emerge in the phase diagram due to the Gaussian 
mode. Thereafter in Section IIVB 41 we conduct a clas- 
sification of the different transitions found in the phase 
diagram. 



A. The Theory 

In this section we develop the perturbative theory to 
derive an effective second order master equation for the 
nuclear system in the vicinity of the Gaussian steady 
state. 

For realistic parameters, the Liouville operator C of 
Eq. (JlJ does not feature an obvious hierarchy, that would 
allow for a perturbative treatment. In order to treat 
the electron-nuclear interaction as a perturbation, we 
first have to separate the macroscopic semiclassical part 
of the nuclear fields. To this end we conduct a self- 
consistent Holstein-Primakoff approximation describing 
nuclear fluctuations around the semiclassical state up to 
second order. 

The (exact) Holstein-Primakoff transformation ex- 
presses the truncation of the collective nuclear spin op- 
erators to a total spin J-subspace in terms of a bosonic 
mode (6 denotes the respective annihilation operator): 



IV. PERTURBATIVE TREATMENT OF THE 
GAUSSIAN MODE 

As seen in the previous section along the segment 
x the system settles in a factorized electronic-nuclear 
state, where the nuclear system can be described as a 
lowering operator eigenstate up to second order in £ = 
J -1 / 2 . Motivated by this result we develop a perturba- 
tive theory based on a self-consistent Holstein-Primakoff 
transformation that enables the description of a class of 
steady states, which generalizes the squeezed coherent 
atomic state solution along x to finite thermal fluctua- 
tions (RSTSS, Appendix [B]) . A solution of this nature 
can be found across the entire phase diagram and we 
show that this treatment becomes exact in the thermo- 
dynamic limit. 

In Section TlVBI we discuss this Gaussian mode across 
the whole phase diagram. Steady state properties of the 
nuclear fluctuations derived from a reduced second or- 
der master equation provide deep insights in the nature 
of the various phases and transitions. Observed effects 



r = y/2J - fctfe b (7) 
I z = tfb - J. 

In the following we introduce a macroscopic displace- 
ment VJ/3 € C (|/3| < 2) on this bosonic mode to account 
for a rotation of the mean polarization of the state, ex- 
pand the operators of Eq. Q and accordingly the Liou- 
ville operator of equation Eq. (jl]) in orders of e = l/\fl. 
The resulting hierarchy in the Liouvillian allows for an 
perturbative treatment of the leading orders and adia- 
batic elimination of the electron degrees of freedom whose 
evolution is governed by the fastest timescale in the sys- 
tem. The displacement /? is self-consistently found by 
demanding first order stability of the solution. The sec- 
ond order of the new effective Liouvillian then provides 
complete information on second order stability, criticality 
and steady state properties in the thermodynamic limit. 

The macroscopic displacement of the nuclear mode 

b^b + y/jp, (8) 
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allows for an expansion of the nuclear operators [Eq. Q] 
in orders of e 



r/J =Vk\jl - e 



f3tf+(3*b n tfb 



<'—{fi + €b) (9) 



where 

J ~ = Vk/3, 



Ji 
J2 



_L[(2fc-|/?|> 
(3*b + (3b\ 



+ 



2y/k 

y/k(3 f r +P*b 2 , >^ 



4 x 



(10) 

(11) 



(12) 



Let P be the projector on the subspace of zero eigen- 
values of Co, i.e., the zeroth order steady states, and 
Q = 1 — P. Since Co features a unique steady state, we 
find Pp = Trs(p) ® p SSl where Trg denotes the trace over 
the electronic subspace and Cop ss — 0. By definition 
it is PCo = CqP = 0. After a generalized Schrieffer- 
Wolff transformation [30], we derive an effective Liou- 
villian within the zeroth order steady state subspace in 
orders of the perturbation 

C cS =ePdP (20) 
+ e 2 (PC 2 P - PCxQC^QCiP) + 0(e 3 ). 

After tracing out the electron degrees of freedom the dy- 
namics of the nuclear fluctuations b are consequently gov- 
erned by the reduced master equation 



&:=Trs(Pp)=Tr s (C oS Pp). 



(21) 



and k = 2 — \(3\ 2 . Analogously, one finds 

2 

i=0 

Jo = m 2 - 1, 
j{ = & + p*b, 

Ji = tfb. 



(13) 

(14) 
(15) 
(16) 



This expansion is meaningful only if the fluctuations 
in the bosonic mode b are smaller than sfj . Under this 
condition, any nuclear state is thus fully determined by 
the state of the bosonic mode b and its displacement j3. 

According to the above expansions Master Eq. ([I]) can 
be written as 



p/J = [C + eCi 



e 2 C 2 



C(e 3 



(17) 



where 



C-oP = 



1 



1 {s- P s+ - -{S+S 



A, 2P : 



p} + ) (18) 
^[5 + (0 + a/2J Q -) + S~{n + a/2J+) 

+ aS + S-J z ,p], 
i[a/2{S + J^ 2 + S-J+ 2 ) (19) 
+ (aS+S- +8u)Jl 2 ,p\. 



The zeroth order superoperator Co acts only on the 
electron degrees of freedom. This separation of timescales 
between electron and nuclear degrees of freedom implies 
that for a given semiclassical nuclear field (defined by 
the displacement /?) the electron settles to a quasi-steady 
state on a timescale shorter than the nuclear dynamics 
and can be eliminated adiabatically on a coarse grained 
timescale. In the following we determine the effective nu- 
clear evolution in the submanifold of the electronic quasi- 
steady states of Co- 



The first order term in e of Eq. (20) can be readily 
calculated 



Tr^PdPp) = -1 [(A) ss b+ (At)„&t,cr] 
where A is an electronic operator 



A =(3*(aS+S- +6lu) 
[(2k- 



(22) 



(23) 



\ 2 )S + 



(p*) 2 S~ 



(A) ss denotes the steady state expectation value accord- 
ing to Cq, which depends on the system parameters 7 
and 17 and on the semiclassical displacement (3 via opti- 
cal Bloch equations derived from Co as described below. 
Eq. (22 1 represents a driving of the nuclear fluctuations 



to leading order in the effective dynamics. Thus for the 
steady state to be stable to first order, we demand 



0. 



(24) 



This equation defines self-consistently the semiclassical 
nuclear displacement /3 in the steady state in dependence 
on the system parameters 7, £1 and Slo. 



The calculation of the second order term of Eq. ( 20 1 is 
more involved and presented in Appendix [Dj We find the 
effective nuclear master equation to second order [M] 



o=2R a { barf - -{tfb,o} 
[tfob-^{bb\o} 



-c ( bob - -{bb.o} 



(25) 



+c* (&W - i{&t&t j0 .}^ 

-i [(/„ + I b + Fpb + (a + B*)b 2 + (a* + B)(tf) 2 ,o] 



with 
B 



aft 



F = 



a 



[(4k- 



2 )(s- 



8Vfc 3 
a((S+S 



(Ak+\ft\ 2 )(ft(S" 
+ 5co/a), 



f3 2 (S + ) i 

-p*(s-) 



and 





poo 

/ df Re ((AUt)A(O)) ) 
Jo 


Ia = 


pOO 

/ dtlm((A\t)A(0)) ss ), 
Jo 


Rb = 


pOO 

/ dtRe((A(t)A\0)) ss ), 
Jo 


h = 


pOO 

/ dtlm((A(t)A^(0)) ss ), 
Jo 


c = 


POO 

/ dt({A(t),A(0)}) ss , 
Jo 


a = 


i f°° dt([A(t),A(0)]) ss . 
Zl Jo 



(26) 
(27) 

(28) 



For a given set of system parameters the coefficients 
defining the nuclear dynamics [Eq. (26 1, Eq. (27 1 and 
Eq. (28)] depend only on the nuclear displacement ft. 



After choosing ft self-consistently to fulfill Eq. (24) in 



order to guarantee first order stability, Eq. ( 25 ) contains 
all information of the nuclear system within the Gaussian 
picture, such as second order stability as well as purity 
and squeezing of the nuclear steady state. Also it approx- 
imates the Liouville operator's low excitation spectrum 
to leading order and thus contains information on criti- 
cality in the system. Eq. (25) therefore forms the basis 



for the subsequent discussion of the RSTSS mode and 
the corresponding phases and transitions in Section 1^ 
In order to calculate the coefficients of Eq. (28). 
we have to determine integrated electronic autocorre- 
lation functions of the type J °° dt (S i (t)S^(0)) 8S and 
f£°dt (S i (Q)S j (t)) ss , where i,j = +,-,z. The dynam- 
ics of single electron operator expectation values are gov- 
erned by the optical Bloch equations derived from Cq 

^(AS)=M(AS), (29) 

where AS := S- (S) ss and S = (5+, S~ , S Z ) T and 



M 



'-(%-ia£$) 


-iO 



.(7 





in* 



(30) 



-2iQ* 

-7 

fi + f Vk(3 and £g is given m 
Eq. (|I4|). The steady state solutions can readily be eval- 

f2*(7 + 2ia£g) 



where we defined f2 



uated 



2i- 



1 2 + AaLf + 8|ft| 2 
I 7 2 + 4a£g 2 
2 7 2 + 4a£g 2 + 8|fi| 2 ' 



(31) 



(32) 



Defining the correlation matrix S = (ASAS^) SS and 
S t = (AStAS^) ss , the Quantum Rc gression Theorem 
[4"T] yields the simple result: 



S\ =(AS AS? 



Se 



(33) 
(34) 



Finally the time integrated autocorrelation functions re- 
duce to the simple expression 



T\ = I dtS t = I ill, 
lo Jo 



dtt 



ML 



S = -M^S, 



dtSl 



(35) 
(36) 



These matrices straightforwardly define the coefficients of 
the effective master equation of the nuclear fluctuations 
Eq. (25). In Appendix D I we provide explicit formulas 



to calculate the relevant coefficients. 



B. Phase Diagram of the Gaussian Mode 

In this Section we use the theory developed above 
to study the RSTSS mode across the phase diagram. 
As outlined in the previous section we first determine 
self-consistently possible semi-classical displacements /?, 
which guarantee first order stability [Eq. (24)]. For 
each of these solutions we determine the effective mas- 



ter equation for the nuclear fluctuations Eq. (25), which 
in the thermodynamic limit contains all information on 
the steady state and the low excitation dynamics and we 
discuss properties like second order stability, criticality as 
well as purity and squeezing of the nuclear steady state. 
Using this information we provide a complete picture of 
the various phases and transitions involving the RSTSS 
solution. 



1. Methods and General Features 

In order to determine the semiclassical displacements 
ft which guarantee first order stability, we show in Ap- 
pendix [C] that Eq. ( 24 ) is equivalent to the semiclas- 
sical steady state conditions. Due to a symmetry in 
the equation, the steady state displacements appear in 
pairs /?+. Any semiclassical displacement (3 can be 
straightforwardly converted to the mean spin polariza- 
tions up to leading order in e according to Eq. (10), 
Eq. (Ill, Eq. pll, and Eq. p2h. In the thermody- 



namic limit the two sets of steady state expectation 
values extracted from /3_ and j3+ share the symmetry 
(±(S x ),(S y ),(S z ),(I x ),±(I y ),±(I z )). In large parts of 
the phase diagram the solution /?_ (f3 + ) displays high nu- 
clear polarization in the same (opposite) direction as the 
the electron spin pumping. We define the corresponding 
quantum states as the normal (anomalous) spin pumping 
mode. 
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The two solutions j3± define two corresponding mas- 
ter equations of the nuclear fluctuations around the re- 
spective semiclassical expectation values according to 
Eq. (25 1. These master equations are subsequently used 



to determine second order stability of the nuclear fluc- 
tuations and, if the dynamics turn out to be stable, the 
steady state properties of the nuclear system. We em- 



phasize that the effective Master Eq. (25) not only can 



be used to determine steady state properties, but also 
reproduces accurately the low excitation spectrum of the 
exact Liouvillian. It thus also describes the system dy- 
namics in the vicinity of the steady state (increasingly 
accurate for large J). 

From Eq. ( 25 ) one readily derives a dynamic equation 



for the first order bosonic moments 



(b) 
(6t) 



((b) 



with 



— (R a — Rt) 

X =I a + h + F, 
£ =a* + B, 



i\ -2i£ 

-(R a - R b ) + ix 



(37) 



(38) 

(39) 
(40) 



where all parameters are functions of the semiclassi- 
cal displacements fj±. This equation of motion - and 
thus the corresponding master equation itself - features 
a fixed point if the eigenvalues of the matrix E have 
negative real part (RefAi^] < 0). Due to the sym- 
metry between /3 + and /3_ one finds that the eigenval- 
ues of the two £ matrices corresponding to j3± fulfill 
Re[Ai,2 (/?+)] = —Re[Ai, 2 (/?-)] such that across the whole 
phase diagram only one solution is stable at a time and 
defines the corresponding phase in the phase diagram. 
Note however, that the unstable solution decays at a rate 
that is second order in e. Preparing the system in this 
state consequently leads to slow dynamics, such that this 
solution exhibits metastability. 

In the following we implicitly choose the stable (3 for 
which the real parts of the eigenvalues of S are negative 
and discard the unstable solution. Fig. [2] displays a selec- 
tion of steady state expectation values in the thermody- 
namic limit across the phase diagram for the stable solu- 
tion. Different expectation values illustrate the different 
nature of phase A and B and show distinct signatures 
of first and second order phase transitions which will be 
discussed in greater detail in Section |IV B 3| and |IVB4| 
The approximate steady state polarizations found in this 
way coincide with the exact values found via diagonal- 
ization techniques to an extraordinary degree (~ 10~ 3 
relative deviation for J=150). Corrections to the pertur- 
bative solutions are of the order 1/J since the first order 
expectation values of the bosonic mode vanish by con- 
struction, since (b) = [Compare Eq. ^ and Eq. (13)]. 
In the thermodynamic limit the perturbative solution be- 
comes exact. 




Figure 2: The system observables of the RSTSS solution in 
the thermodynamic limit show clear signatures of first and 
second order transitions, (a) The nuclear polarization in z- 
direction {I Z /J) 3S switches abruptly from minus to plus at the 
phase boundary b. (b) The electron polarization in z-direction 
{S x ) ss shows a similar discontinuous behavior along b. (c) The 
nuclear polarization in ^-direction changes smoothly across 
the phase boundary b. Along the segment x (a; = Uo, fi < fio) 
the nuclear field in x-direction builds up linearly to cancel the 
external driving, (d) The electron polarization in z direction 
also does not show signatures of the first order transition b. 
Along segment x the electron is fully polarized in -z direction 
up to the second order critical point (uo, fio), where it changes 
non-analytically (see also Fig. j6j). 



The two eigenvalues of E are typically of the form 
X12 = a±ib (except in region D which will be discussed 
below) and define the complex energy of the mode. In 
this case the matrix E contains all information on the low 
excitation spectrum of the Liouvillian which is approxi- 
mated by multiples of the mode energies within the per- 
turbative treatment |55j . The low excitation spectrum 
contains information about criticality of the system and 
the dynamics in the vicinity of the steady state, and will 
be used to discuss and classify the different transitions 
in the phase diagram. In particular the eigenvalue of E 
with largest real part approximates the ADR in the ther- 
modynamic limit in those regions of the phase diagram 
where the Gaussian mode is responsible for the lowest ex- 
citations in the Liouvillian spectrum (only in the region 
of bistability C this is not the case). 

The ADR according to the perturbative descriptions 
based on Gaussian modes is displayed in Fig. [3] It is used 
to study the transitions involving the Gaussian mode in 
the thermodynamic limit. The ADR vanishes along a 
line b indicating a phase boundary separating the normal 
and anomalous spin pumping phase, which is described 
in Section |IV B 4| Furthermore a non-analyticity of the 
ADR at a finite value defines region D, which character- 
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izes a dynamical phase and is explained in Section IV B 3 
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Figure 3: Asymptotic decay rate (ADR, cf. text) for 7 = a 
within the perturbative framework. Along b the ADR van- 
ishes non-analytically indicating the stabilizing and destabi- 
lizing of the modes of region A and B, respectively. 6 is a first 
order phase boundary culminating in a second order critical 
point at (a;o,fio)- From here region D opens which is char- 
acterized by a non-analyticity in the ADR at a finite value. 
This indicates a change in the dynamic properties of the sys- 
tem which can not be detected in steady state observables. 
Within D the system shows an over damped behavior in the 
vicinity of the steady state. 

The dynamical matrix of the first order moments £ 
provides information on the stability of the semiclassi- 
cal solutions, the criticality of the Liouvillian and the 
non-analyticities of region D. In order to understand the 
character of the solutions in the different regions of the 
phase diagram we consider next the steady state covari- 
ance matrix (CM) of the bosonic system. For a quadratic 
evolution like the one of Eq. ( 25 ) the steady state covari- 
ance matrix contains all information on the state. We 
deduce the effective temperature and the squeezing of 
the nuclear spin system, which connects to criticality in 
the system. 

For a one-mode system with vanishing displacements 
(x) and (p) [in the steady state of Eq. ( 25 1 this is always 
the easel the CM is defined as 



r = 



2{x 2 ) 
2 (pa;) -I 



2(xp) - 
2(P 2 ) 



(41) 

with the usual definitions x = + and p = 

-^Xb — tf). Using Eq. (25) we straightforwardly cal- 
culate the steady state covariance matrix r ss across the 
phase diagram. As T = T T > 0, T is symplectically 
diagonalizable, with 



r = do 



M 2 
M- 2 



(T 



(42) 



where O is orthogonal with det(O) = 1. For a single 
mode, D > 1 and M > 1 are real numbers. While D is 
a measure of the purity of the state \Tr{p 2 ) = l/y/\T\ = 
l/£>], the smallest eigenvalue of T, X m [ n = DM~ 2 de- 
termines the amount of squeezing in the system [42] . 
Amin < 1 indicates squeezing in the bosonic mode. For 
M = 1, the covariance matrix Eq. (42 1 describes a ther- 



mal state of the bosonic mode and D can be straightfor- 
wardly associated to a dimensionless effective tempera- 
ture 



T e ft = In 



Vd- 1 



1 



(43) 



This definition is also meaningful for M > 1, since the 
squeezing operation is entropy-conserving. T e g is also a 
measure for the entropy of the spin system, as to leading 
order it is connected to the bosonic mode via an unitary 
(i.e., entropy-conserving) transformation. The effective 
temperature of the different phases will be discussed be- 
low in Sections |IVB2| & |IVB3| [cf. Fig. [F] . 



We stress the point that all properties of the covariance 
matrix derived within the second order of the perturba- 
tive approach are independent of the system size J. In 
particular, the amount of fluctuations (i.e., the purity) in 
the state does not depend on the particle number. In or- 
der to self-consistently justify the perturbative approach, 
D has to be small with regard to J. This implies that in 
the thermodynamic limit J —¥ 00 the perturbative results 
to second (i.e., leading) order become exact. The inverse 
purity D is displayed in Fig. [4] a) . Except for for a small 
region around the Gaussian phase boundary b the fluc- 
tuations are much smaller than J = 150, which justifies 
the validity of the perturbative approach and explains 
the excellent agreement with the exact diagonalization 
for this system size. 

The squeezing A m ; n in the auxiliary bosonic mode does 
not necessarily correspond to spin squeezing in the nu- 
clear system. In order to deduce the spin squeezing in the 
nuclear system from the squeezing of the bosonic mode a 
transformation according to Eq. ( 11 ) and Eq. (151 is nec- 



essary. In Appendix A 1 we show that for \/3\ < 1 Eq. ( 11 ) 



can be reformulated to connect the spin fluctuations to 
a squeezed and rescaled bosonic mode 



Jf = VW^W)S\r)bS(r) 



(44) 



where S(r) — e ( - r b rb+ ^ 2 is the squeezing operator and 
cosh(r) =fi = (2fc-|/3| 2 )/[2^2fc(l - |/3| 2 )] and sinh(r) = 

-^ = /3 2 /[2 v /2fc(l-|/3| 2 )]. 

Thus squeezing A m j n of the mode b does in general not 
imply reduced spin fluctuations in a direction orthogonal 
to the mean spin polarization since the transformation 
between spin fluctuations and b involves a squeezing op- 
eration itself and a scaling by a factor < \/2(l — |/3| 2 ) < 

V2. 

In general we thus have to apply a more involved 
squeezing criterion. In [35] it was shown that for sys- 
tems of N spin- 1/2 particles and for all directions n the 
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Figure 4: Properties of the steady state covariance matrix V S3 
[Eq. ( 42 1] . a) The fluctuations D are low in most parts of the 
phase diagram except for a small wedge around the Gaussian 
phase boundary, b) Fluctuations D along the line u) = 1.5 Uo 
[green line of a)]. The phase boundaries separate a mode 
with low fluctuations (enlarged in the inset), from a mode 
with large fluctuations. For large Q fluctuations increase, and 
the system eventually approaches a fully mixed state, c) The 
squeezing measure C (c.f. text) in the thermodynamic limit. 
C approaches 1 at {uiq, ^o) indicating diverging entanglement 
in the system, d) C along the line to = loo (solid line). The 
red circles indicate the the squeezing parameter 1 — ff^ = 

1 - sjl - (ftM,) 2 (cf. text). 



quantity 



1 



1 

J2 



In) 2 < 1, 



(45) 



signals entanglement if > for some direction n. 
Moreover, (A/S) < J/2 indicates a generalized spin- 
squeezing of the state [55] . 

In the following we will use the quantity C = 
max{0, Cft \ri G R 3 } to investigate squeezing and bipar- 
tite entanglement in the nuclear system. In order to cal- 
culate Cn we reconstruct the approximate nuclear opera- 
tors according to Eq. ^ and Eq. ( 14 ) from the semi- 
classical displacement /3 and evaluate the expectation 
values according to the steady state covariance matrix 
Eq. (41 ). Finally we maximize Cn with regard to all pos- 
sible directions n to obtain C . The results arc discussed 
in Section HV B 41 As discussed in more detail in the next 
Section, the fact that C — > 1 as ft — > f^o on the line seg- 
ment x indicates a diverging entanglement length in the 
sense that 0(1/(1 — C)) = 0(\/j)-particle entanglement 
is present [37]. 



2. A Second Order Phase Transition: The Segment x 

The segment x at uj = uj (Fig. [lj represents a very 
peculiar region in the phase diagram, where the solution 



below the critical point can be constructed analytically as 
seen in Section flH B The electron and nuclear system de- 
couple, resulting in a zero entropy product steady state. 
A nuclear polarization builds up to cancel the external 
driving up to the point of maximal Overhauser field (fio)- 
At this point squeezing and entanglement in the system 
diverge, indicating a second order phase transition. In 
the following we exemplarily employ the formalism de- 
veloped above along this line to obtain further insight 
about the criticality at (u!o,flo). We calculate the ana- 
lytical steady state solution as well as the effective master 
equation governing the nuclear fluctuation dynamics in 
its vicinity. We find that here the spectrum of the Liou- 
villian becomes continuous (implying a closing gap) and 
real. At the same time the creation operators of the ele- 
mentary excitations from the steady state turn hermitian 
giving rise to diverging spin entanglement. 
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Figure 5: The ADR (7 = a) for J = 50, 100, 150 (broken 
lines) in comparison with the perturbatively calculated (solid 
line, cf. Section |IVB"2"| ) along uj = cjo- For finite systems one 
finds an avoided crossing at Qo- The size of the gap reduces 
with the system size until it closes in the thermodynamic limit 
(solid line). Below f^o the ADR in the thermodynamic limit 
is given by Eq. (52 1. 



The first order stability condition Eq. (24 1 is fulfilled, 
if Q = [compare Eq. fell) and Eq. fef], which yields 
the possible semiclassical steady state displacements 



HcP = -fi/Qo 



(46) 



<= 



/3± = -/liv/Hiw, 



corresponding to a normal ('— ') and anomalous ('+') spin 
pumping mode, respectively. 

Next, we explicitly calculate the second order correc- 
tive dynamics of the nuclear degrees of freedom for the 
normal mode. The vanishing of the effective driving 
= forces the electron in its dark state - implying 
(S + ) ss = (S-) ss = (S + S-) ss = - and directly yields 
B = F = [Eq. ([26]) and Eq. p7b]. The remaining 
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constants can be calculated as described above and in- 
troducing new bosonic operators (for the normal mode 
P = P-<1) 



d = fib + v$ , 



(47) 



with 



2k 



2 v /2fc(T~ 

/? 2 



2V2fc(l -|/?| 2 )' 



(48b) 



one finds the effective evolution of the nuclear fluctua- 
tions given as 



<j=r, 



cfr 



dad) - ~{rfd,a}j 



(49) 



i [O c gd^d,a] 



with 



c ff = a Im 



7 + i 2a(|/3| 2 -l) 



(1-|/3| 2 ). (51) 



d and fulfill boson commutation relations, since 



Eq. (47) defines a symplectic transformation 



1). The eigenvalues of the dynamical matrix £ associ- 
ated to Eq. (49) are straightforwardly given as Ai,2 = 



— r c ff/2 ± iOcft ■ The real part - representing the ADR of 
the system in thermodynamic limit (compare Fig. [5]) - 
is always negative, indicating the stability of the normal 
spin pumping mode (/?-). In an analogous calculation 
one shows that the semiclassical solution /3+ > 1 is not 
stable to second order since the eigenvalues of S have a 
positive real part, i.e. the fluctuations diverge, violating 
the initial assumptions that the mode b has to be lowly 
occupied. 

Selected steady state expectation values derived from 
the stable displacement f3_ to leading order in J (i.e., 
in the thermodynamic limit) are displayed in Fig. [6] 
Already for J = 150 we find excellent agreement be- 
tween the perturbative and exact mean polarizations. 
The semiclassical nuclear field builds up to exactly can- 
cel the external magnetic field forcing the electron in 
its dark state \\.) along x and thus realizing the model 
of cooperative resonance fluorescence [7] even for weak 
dissipation 7 < a [compare Eq. ([5|]. This solution is 
only available if < fio (defining segment x), i.e., up to 
the point where the nuclear field reaches its maximum. 
At this point the system enters a new phase of anoma- 
lous spin pumping (described below) and the steady state 
properties change abruptly. 

Inserting solution /3_ in the coefficients of Master 
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Figure 6: Electron inversion (S z ) and the nuclear field in 
x direction (I x ) along cj = ojo, in the thermodynamic limit 
according to the perturbative theory (circles) in comparison 
with the numeric values from exact diagonalization for a fi- 
nite system of J = 150 (solid lines). The perturbative theory 
shows excellent agreement with the numerics. Further the 
numerically determined electron inversion and the expecta- 
tion value of the inhomogeneous nuclear operator (A x ) are 
displayed for a model of two inhomogeneously coupled nu- 
clear shells (gi = 2 172) of size Ji,2 = 8 (dashed lines) and for 
5 inhomogeneously coupled nuclear spins (dotted lines) are 
displayed (discussion see Section VI I. 



Eq. (|49j) yields 
r off = 2a 2 Re ^ 

©eff = a 2 Im 



7- aa^/\ - {n/n Q y 



7-i2ay / l - (^M)) 2 



(52) 

y/l - (f>/U) 2 . 

(53) 



In the close vicinity below the critical point tto the real 
part of the gap in the Liouvillian's spectrum closes as 



r cff «2- v /i-(^M) s 

7 



and the imaginary part as 



|e cff |«2— [i-(fi/n ) 



(54) 



(55) 



indicating criticality. Fig. [5] displays the ADR along 
u) = luq in the thermodynamic limit (which is given on 
the segment x by Eq. ( |52[ )) and for finite systems. It 
displays an avoided crossing at fio with a gap that van- 
ishes in the thermodynamic limit. This closing of the gap 
coincides with diverging timescales in the system, which 
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renders the model more susceptible to potential perturb- 
ing effects, a phenomenon well known in the context of 
criticality 39 . 



In contrast to the general form Eq. (251, Eq. (49) con- 



tains only one Lindblad term and the dynamics drive 
the system into the vacuum |0<j) of the squeezed mode 
d. As the system approaches the critical value Q = 
(i.e., /3_ = —1) the mode d adopts more and more a 
p = -^{b — b') like character and thus the squeez- 
ing of this mode's vacuum increases. The (in general 
complicated) transformation between the squeezing of 
the bosonic mode b and the spin operators (cf. Sec- 
tion [IVBT) can readily be established along x, since the 
operator d is trivially related to the spin operators [cf. 
Eq. pj] 



1 



[(2k 



2 )6-/3 2 6t] 



= v/2(l- W)W> + vtf) 

= V2(i-|/W 



(56) 



The fluctuations in y-direction, for example, are conse- 
quently given as 



where pd = ^^(d — <v). One readily shows that 



(AI 2 ) = J{jf) = J{l-WW d ), 



(57) 



(58) 



up to order 0(1) and we used (d) = in the steady 
state. In the p-vacuum |0 P ) it is (p 2 ) = 1/2, such that 
we evaluate 



tt y =2{Ml)/\{T)\ 

= 2(l-\P\ 2 )(Pl)=Jl- 



(59) 



where we used \(I)\ = J and inserted the semiclassical 
displacement /3_. 

This is the same result we derived in Section UlI Bl and 
Appendix |A 1| by constructing approximate eigenstates 
of the lowering operator J - and along x we find that 
C sa 1 — , as shown in Fig. jij d) . Note that here e y 

is orthogonal to the direction of the mean spin (7) . This 
allows us to deduce that 0(y/J) nuclear spins must be 
entangled close to the critical point, which establishes a 
'diverging entanglement length' in this system. To see 
this, we employ a variant of the criterion Eq. (|6| as dis- 
cussed in [36]. There, it was shown that £| < 1/k sets 
a lower bound of N£~ 2 on the quantum Fisher informa- 
tion Fq of the state. In [37 it was shown that for states 
containing at most fc-particle entanglement, Fq is upper 
bounded by Nk. Consequently, the values of £f obtained 



close to the critical point (cf. Eq. ( 59 ) and Appendix A 1 ) 
imply that at least 0(vJ)-particle entanglement must be 



present. Note that the bosonic description does not al- 
low to describe the range £§ = 0(1/ J), i.e., k = O(J), 
where the fluctuations become larger than the expansion 
parameter. 

The nuclear squeezing and entanglement in the system 
diverges approaching the critical point, as the Lindblad 
operator d (defining the steady state |0d)) becomes more 
and more p-like. The fluctuations in y-direction tend to 
zero, while at the same time - due to the Heisenberg un- 
certainty relation - the steady state is in a superposition 
of an increasing number of I z eigenstates. Since in a sys- 
tem with infinite range interactions (as the one we are 
considering) there is no obvious definition of a coherence 
length, the range of the involved I z eigenstates can be 
considered as an analogous concept. 

At the critical value O = the symplectic transfor- 



mation Eqs. (47) becomes ill defined (d becomes a p-like 
operator) while both the dissipation rate and the mode 
energy tend to zero. While the coefficients in Eqs. (48) 



diverge, the total master equation is well defined [due to 
the factors (1 — \/3\ 2 ) in r o g] and straightforwardly can 
be written as 



. a 2 / 



(60) 



The Liouville operator's spectrum is real and continu- 
ous with hermitian creation operators of the elementary 
excitations. 

We stress the point that along segment x in the phase 
diagram highly dissipative dynamics drive the system in 
a pure and separable steady state with zero effective tem- 
perature T c ff — [cf. Fig. [7]b)]. At the critical point flo 
the steady state changes its nature abruptly as the sys- 
tem enters a high temperature phase. 

Furthermore we remark that this steady state has no 
relation to the system's ground state. This is in con- 
trast to the extensively studied Dicke phase transition 
P21 |43l |44] where the steady state is in close relation 
to the Hamiltonian's ground state (in fact in the normal 
phase it is identical). In the present model dissipation 
drives the system to a highly excited state of the Hamil- 
tonian and the observed critical phenomena are discon- 
nected from the Hamiltonian's low excitation spectrum. 

We have seen that at the critical point (u>q,D,q) the 
gap of the Liouville operator's spectrum (in both real 
and imaginary part) closes in the thermodynamic limit 
[Eq. (54) and Eq. (55)]. Approaching the critical point 



the steady state fluctuations become more and more 
squeezed due to the increasing j5-like character of the 
mode d. The spin squeezing close to the critical point 
[Eq. ( 59 1] can be interpreted as a diverging coherence 



length in a system with infinite range interactions (the 
electron mediates interactions between remote spins). 
These are clear indications for a second order phase tran- 



sition, which will be formalized in Section IV B 4 
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3. Phases 

In the present Section we study the different phases 
of the system, which involve the RSTSS solution (A, B 
and D) using the analytic tools developed above. By 
construction, the RSTSS solution describes steady states 
where the electron and nuclear state factorize to lead- 
ing order in the system size and the nuclear system is 
found in a fully polarized and rotated state with Gaus- 
sian fluctuations, which are fully characterized by their 
effective temperature and squeezing. Fig. [2] displays dif- 
ferent steady state observables of the Gaussian solution 
determined via the formalism described above to leading 
order in the thermodynamic limit. 



RSTSS now is in a high-temperature state. For larger 
electron driving the temperature increases until eventu- 
ally the Gaussian description breaks down (as D oc J) 
and for f2 — > oo the system is found in a completely 
mixed state [compare Fig.[4]b)]. 

In the upper half of the phase diagram (lj > u>o) phase 
A changes abruptly into phase B at the boundary b and 
certain steady state spin observables [(I z ), (S x ) [Fig.[2]a) 
& b)] and (I y ) (not displayed)] show distinct features of a 
first order phase transition, changing sign as the normal 
(anomalous) mode destabilizes (stabilizes). This transi- 
tion is discussed in greater detail in the following Sec- 
tion IV B 4 Following this boundary towards the criti- 
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Figure 7: Effective temperatur T e s of the gaussian mode. 
Temperatures T c s > 6 are cut off, as the temperature di- 
verges along the phase boundary b. a) The first order phase 
boundary b separates the low temperature phase A from the 
high temperature phase B. b) T g along ui — ujo- On segment 
x the system is in a zero entropy state (T G ff=0). Above the 
second order critical point £1 > S7o the system enters a high 
temperature phase. Here the temperature rises with increas- 
ing driving strength. 

In phase A the system is characterized by normal spin 
pumping behavior. Only the semiclassical displacement 
/?_ (normal mode) leads to a dynamical matrix E that 



cal point (wo,^o) the two phases become progressively 
more similar. Below the critical point (w < ujq) there 
is no clear distinction between the normal and anoma- 
lous spin pumping mode anymore, a phenomenon known 
from thermodynamics as supercriticality. Phase A trans- 
forms continuously to phase B in this region. Close to the 
critical point, supercritical media typically respond very 
sensitive to the external control parameters of the phase 
diagram (e.g., temperature or pressure) [38]. In our sys- 
tem we observe that small changes in the parameter co 
leads to large changes in electron spin observables. 

Next, we consider the third region associated with the 
RSTSS solution, region D. We will find that this region 
differs from the previous ones by the fact that it cannot 
be detected in the system's steady state but rather in 
dynamical observables. 

The eigenvalues of the dynamic matrix E can be calcu- 
lated as Ai,2 = — (R a ~ Rb) ± 2a/4|£| 2 — x 2 and provide 
information on the approximate low excitation spectrum 
of the Liouvillian. We can distinguish two cases for the 
low excitation spectrum, which differ only in the Hamil- 



tonian properties of Eq. (251 (fully determined by x an d 
£ [Eq. p9l) & Eq. pi]). In the first case the quadratic 



bosonic Hamiltonian can be symplectically transformed 
to be diagonal in a Fock basis (i.e., of the form oc b^b). 
has negative reIlparts~of its e7genvalue^rwMe for /3+ the This is the case if ^ > As a consequence the two 



eigenvalues have positive real parts, indicating the insta- 
bility of that mode in second order. The nuclear system 
in the normal mode settles in a state highly polarized in 
-z direction following the direction of the electron spin 
pumping [Fig. [2ja)]. Meanwhile, increasing the external 
driving f2 and approaching the phase boundary b, a nu- 
clear field in x direction builds up, but only along x it 
can fully cancel the external driving [Fig. |2^c)]. There- 
fore, in general the electron spin aligns more and more 
with the external field [Fig. [2jb,d)]. Furthermore, the ef- 
fective temperature (and thus the entropy) of the phase 
is low, as displayed in Fig. [7] a) . 

In region B in contrast, f3 + is the only stable solution, 
defining the phase of anomalous spin pumping behav- 
ior. The nuclear system now shows strong population 
inversion, i.e., the nuclear polarization is in direction op- 
posite to the external pumping (z). In the same way 
the electron now aligns in opposite direction to the ex- 
ternal driving field (x). Also, in contrast to phase A, the 



eigenvalues of E have an identical real part and imaginary 
parts ±2\J 'x 2 — 4|£| 2 . In the second case the Hamiltonian 
transforms symplectically into a squeezing Hamiltonian 
cx (6 t2 + b 2 ). Here one finds x 2 < 4|£| 2 , such that the 
eigenvalues become real and symmetrically distributed 
around —(R a — Rb)- In region D in Fig. [l] we find the 
effective Hamiltonian for the nuclear fluctuations to be 
symplectically equivalent to a squeezing Hamiltonian. 

Fig. [8] shows the ADR exemplarily along the line w = 
0.5 w (@ in Fig.[j} calculated according to the pertur- 
bative theory and via exact diagonalization, respectively. 
The perturbative theory approximates accurately the low 
excitation spectrum of the Liouvillian. We find that in 
region D the ADR splits up, when the coherent part of 
Eq. ( 25 ) changes to a squeezing Hamiltonian. As men- 



tioned above this non-analyticity occurs at a non-zero 
value of the ADR and thus does not leave signatures in 
the steady state behavior. The steady state transforms 
smoothly along (n). However the nature of dynamical 
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Figure 8: The ADR and the imaginary part of the respective 
eigenvalue (7 = a) for J = 150 (solid lines) in comparison with 
the perturbatively calculated value (dots) along © of Fig. [3] 
In the region where the coherent part of Eq. ( |25[ ) is a squeezing 
Hamiltonian, the ADR (i.e. real part of the lowest Liouvillian 
eigenvalue pair) splits. At the same time the imaginary part 
of the lowest eigenvalue pair vanishes (black dashed lines), 
indicating that the system is overdamped. 



observables change within region D as the system dis- 
plays anomalous behavior approaching the steady state. 
The splitting of the ADR coincides with the vanishing 
of the imaginary part of the lowest non-zero Liouvillian 
eigenvalues. Thus the system is overdamped in D. Per- 
turbing the system from its steady state will not lead 
to a damped oscillatory behavior, but to an exponential, 
oscillation-free return to the steady state. 

The blue area in vicinity to region D in Fig. [3] does not 
represents a new phase but is another interesting feature 
of the system. Here, the ADR exceeds the value at ft = 
by a factor ~ 3. For £7 — the model describes the 
standard spin pumping setting. Large gaps in the low 
excitation spectrum indicate the possibility to improve 
the effective spin pumping rate (remember that also in 
this region the steady state is fully polarized, however 
not in — z direction as it is the case for the normal spin 
pumping configuration £7 = 0). Indeed simulations show 
that starting from a fully mixed state, the system reaches 
the steady state faster than in the standard setting (£2 = 
0). This feature becomes more distinct in systems, where 
the electron pumping rate 7 is limited. For 7 = 0.1a the 
time to reach the fully polarized steady state from a fully 
mixed state is shortened by a factor ~ 6. 



4- Transitions 

In this Section we consider the transitions involving 
the RSTSS solution in greater detail providing a classifi- 
cation in analogy to quantum phase transitions in closed 
systems (compare Section |ll|) . 



As seen in the previous Section, certain steady state 
observables show clear signatures of a first order phase 
transition at b (Fig. [2]). In order to understand this sharp 
transition we consider the ADR exemplarily along path 
(T) in Fig. [9J The broken lines represent numeric results 




Figure 9: The ADR (7 = a) for J = 50, 100, 150 (broken 
lines) in comparison with the perturbatively calculated (solid 
line) along of Fig. [3] The vertical black lines indicate 
the asymptotic boundaries of the region of bistability. In the 
whole region the ADR tends to zero in the thermodynamic 
limit due to the appearance of a non-gaussian stable mode. 
Inset: The next higher excitations in the spectrum for J = 
150 display equidistant splittings in regions far from the region 
of bistability. This is an indication for the bosonic character 
of the steady state, which is exploited in the perturbative 
approach. 

of exact diagonalization of the Liouvillian for J = 50, 100 
and 150, while the solid line indicates the result of the 
perturbative approach. As described in Section |TVB l| wc 
implicitly choose the semiclassical displacement /?_ (for 
il < 1.5Q ) or /?+ (for Q > 1.5O ) for which the ADR 
is negative, indicating a stable solution. For increasing 
system size the ADR is increasingly well approximated 
by the perturbative solution. 

We stress the point that the red line represents the 
first Gaussian excitation energy only. However, within 
the region of bistability (indicated by two vertical bars 
and discussed below in Section [V|, a non-Gaussian mode 
(discussed below) is responsible for additional excitations 
in the exact spectrum. The Gaussian mode eigenvalue 
(red line) in this region is reproduced approximately by 
higher excitations of the exact spectrum (not displayed) . 
The perturbative theory is still correct within the region 
of bistability but, as expected, it misses all non-Gaussian 
eigenstates of the exact Liouvillian. 

At the boundary 6 ( (2 « 1.5 flo ) the g a P m the 
real part of the spectrum of the Liouvillian closes non- 
analytically, indicating critical behavior. This observa- 
tion is supported by the effective temperature (and thus 
the fluctuations in the system) , which is increased in the 
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vicinity of the boundary b, and diverges at the boundary 
[Fig. [7] a) & Fig. [I] a)]. The vanishing of the ADR at b 
(i.e., the vanishing due to the RSTSS solution) can be 
observed at finite J (dashed lines in Fig. [9]) and is not a 
feature appearing in the thermodynamic limit only. The 
position of this closing of the gap - which in the ther- 
modynamic limit (solid line) is found at Q w 1.5 £Iq - is 
shifted for finite system sizes to lower drivings f2. 

The origin of this closing of the Liouvillian gap be- 
comes more transparent if we take the mode energy of 
the respective metastable solution into account. 



In Fig. 10 a) the complex energy of both the stable and 
unstable mode are displayed (i.e., the first eigenvalue of 
the matrix £ [Eq. ([37])] ) . The normal spin pumping mode 




Figure 10: Complex energy of the two modes corresponding 
to the semiclassical solutions /3± for 7 = a. The solid line in 
the non-shaded area represents the ADR of Fig. [9] and Fig. [5] 
respectively, a) Along (ui = 1.5 uiq). The eigenvalues 
miss each other in the complex plane. The real parts cross 
directly, b) to = uiq. The eigenvalues degenerate asymptoti- 
cally (in both real and imaginary parts) at the critical point. 
This closing of the gap originates from an avoided crossing in 
finite systems with the relevant gap vanishing in the thermo- 
dynamic limit (see also Fig. [5| 

(/?_; blue lines) is stable [ReA(/3_) < 0] up to the criti- 
cal point where it destabilizes and the anomalous mode 
appears (/3+; red lines). At the critical point the two so- 
lutions are macroscopically different /3_ 7^ /3 + and their 
energy [Im(A/3±)] is distinct across the transition [dotted 
lines in Fig. 10 a)]. Although the projection of the eigen- 



values on the real axis vanishes at the critical point for 
both modes (indicating the stabilizing / destabilizing of 
the modes) the eigenvalues pass each other in the com- 
plex plane at large distance. There is no degeneracy in 
the spectrum of the Liouvillian at the critical point and 
consequently there can be no mixing of the two modes; 
the real parts of the eigenvalues cross directly without 
influencing each other. Except for the change in stabil- 
ity the modes do not change their character approaching 
the phase boundary and no diverging correlations (in- 
dicated by the squeezing parameter C) can be observed. 
Together with the discontinuous change in system observ- 
ables such as mean polarizations we classify this Gaussian 
transition as of first order. 

Second, we consider the transition along uj = ujq (in- 



cluding the line segment x). In contrast to the situation 
before we find that the semiclassical displacements /3+ 
and /3_ merge approaching the critical point such that 
the two modes become asymptotically identical at £Iq 
[Eq. ( 46 )] . Approaching the critical point, the eigenvalues 
of the two modes tend to zero (both the real and imag- 
inary parts), causing the gap of the Liouvillian's spec- 
trum to close [Fig. 10 b), Eq. ( [54] ), Eq. (J55"j)]. As we have 
seen in Section IV B 2 at (wo,^oj the spectrum becomes 
real and continuous signaling criticality. The perturba- 
tive treatment intrinsically is a description in the ther- 
modynamic limit. If wc consider the exact spectrum we 
indeed find an avoided crossing due to the mode mixing 
at the critical point with a gap that is clos ing for J — > 00 
(cf. Fig. [5]). As we discussed in Section IV B 2 the el- 



ementary excitations become p-like, causing a diverging 
coherence length in the system [indicated by the diverg- 
ing squeezing parameter C in Fig.[4]jc,d)]. Together with 
the continuous but non-analytical change of the mean 
polarizations these properties classify the point (J7o,wo) 
as a second order transition. 



V. REGION OF BISTABILITY: 
NON-GAUSSIAN SOLUTION 



As noted in Section III B| along the Gaussian boundary 
b extends a region of bistability [C in (Fig. [I )] - culmi- 
nating in the critical point (T1 0i ujq) - in which a second 
stable solution appears. Within the perturbative frame- 



work from Section IV this highly non-Gaussian solution 
could not be detected because it features large fluctua- 
tions of the order of the system size J. In the following 
we use numerical techniques to construct and study this 
mode for finite systems. In the thermodynamic limit the 
ADR tends to zero within C, such that there exists a two 
dimensional subspace of steady states. Here we find two 
independent, physical solutions within the two dimen- 
sional kernel of the Liouvillian, one of which will turn out 
to be the Gaussian normal spin pumping mode described 
in Section [IV} We analyze the nature and properties of 
the other, non-Gaussian solution, exemplarily along the 

line uj = 1.5 u> (© in Fig. [j]). 

Fig. [9] displays the ADR for different particle numbers. 
Within the indicated region of bistability (the black ver- 
tical lines represent the boundaries c and b, respectively) 
the ADR tends to zero with increasing particle number. 
Already for J = 150 one finds a small region, where the 
ADR is small enough (of the order of 10 _6 a) that one can 
construct two linearly independent (quasi) steady state 
solutions. Although we find the eigenmatrix p\ associ- 
ated with the ADR to be non-positive and traceless (the 
latter being a consequence of C being the generator of 
a trace preserving map) we can linearly combine it with 
the true steady state po to obtain two linear independent, 
positive solutions with trace one, p\ (corresponding to 
the normal spin pumping mode) and p up . These solu- 
tions span the two dimensional space of steady states in 



17 



that region. 




Figure 11: Diagonal elements p(m) — (mj p\m) of the nu- 
clear density matrix in z- basis (I z \ m) = m\m)) across the 
region of bistability for J = 150, 7 = a. In the bistable 
regions two quasi stable modes - the Gaussian normal spin 
pumping mode (lower branch; p\ ) and a non-Gaussian (up- 
per branch; p up ) - coexist. The blue dots (red diamonds) in 
the plane indicate the average polarization in z-direction (I z ) 
for the lower (upper) solution. 



Fig. 
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illustrates the solutions p\ and p up around the 
bistable region in an equally weighted mixture. The den- 
sity matrices are represented by their diagonal elements 
in the I z basis. In the plane the blue dots (red diamonds) 
represent the polarization in z-direction (I z ) of the lower 
(upper) solution p\ a (p U p)- Coming from below the crit- 
ical region (£1 < 1.15 VLq) the nuclear system is found 
in the Gaussian normal spin pumping mode, fully po- 
larized, slightly rotated away from the — z direction and 
with fluctuations of the order of \fj . This Gaussian solu- 
tion persists within the critical region where it becomes 
noisier until eventually - approaching the right boundary 
b at Q, = 1.5 f2o ~ it destabilizes. In the thermodynamic 
limit the lower solution is stable up to the right boundary, 
where a first order transition occurs and the anomalous 
spin pumping mode appears. Approaching boundary b 
from above (ft > 1.5 f^o) this mode transforms into a 
non-Gaussian solution, which - in contrast to the coex- 
isting normal mode - features fluctuations of the order of 
J, is not fully polarized and shows large electron-nuclear 
and nuclear-nuclear connected correlations. Approaching 
the left boundary c at O = 1.15 fio this mode destabi- 
lizes eventually as the ADR becomes finite again and the 
normal mode is the only stable solution in the system. 

The bistable behavior of the system in region C bears 
close resemblance to the phenomenon of optical bista- 
bility for saturable absorbers [45] . where connections to 
phase transitions have been established [35]. In this re- 
gion the system displays strong hysteretic behavior. Re- 
cent experiments in quantum dots, realizing a setting 
close to our model system display distinct signatures of 
hysteresis upon application of an external driving field 



on the electronic spin [M] [55] . Our results suggest the 
observed optical bistability in central spin systems as a 
possible pathway to understand these experimental re- 
sults, which will be a subject of further studies. 



VI. IMPLEMENTATIONS AND EXTENSIONS 
OF THE MODEL 

In the present section we discuss potential physical re- 
alizations of the Master Eq. ([I]) and address certain as- 
pects of an extension of the model for inhomogeneous 
hyperfine couplings. 

As mentioned above the model we study is a generic 
central spin model with various potential physical imple- 
mentations. The most prominent ones represent singly 
charged semiconductor quantum dots, where the elec- 
tron spin couples to the nuclear spins of the host mate- 
rial |23L I34| , and diamond nitrogen vacancy centers cou- 
pled to either nuclear ( 13 C spins of the host material) or 
electron (e.g., nearby nitrogen impurities) spin ensembles 
[461 147] . Recently diamond nano-crystals containing sin- 
gle NV centers coated with organic molecule spin labels, 
which are dipole coupled to the NV center spin have been 
manufactured }48| . 

NV centers represent a natural realization of the Mas- 
ter Eq. ([I]) . Their ground state consists of three spin sub- 
levels (of spin projection quantum number m = 0,±1) 
featuring a zero field splitting due to anisotropic crystal 
fields of 2.88 GHz 46 1. In a static magnetic field this 
zero field splitting can be compensated for and one of 
the transitions (e.g., m = O 1) is brought into near 
hyperfine resonance with the ancilla spin system, defin- 
ing an effective two-level system. Since the m = level 
does not carry a magnetic moment, the hyperfine inter- 
action of the effective two level system and the ancilla 
system takes the anisotropic form of Eq. Q. Potential 
countcrrotating terms of the dipole-dipole interaction are 
neglected in the static magnetic field in a rotating wave 
approximation. Optical pumping of the electron spin in 
the m = spin state and resonant driving (either by op- 
tical Raman transitions or radio frequency fields) realizes 
Master Eq. Q [3D]. 

In general the hyperfine interaction in such a setting 
will not be homogeneous and the truncation to a sym- 
metric subspace of total spin J is not justified. In the 
following we consider an extension of the model taking 
into account the inhomogeneous nature of the hyperfine 
coupling in a shell model. Along x we show that up to 
the critical point steady states can be constructed an- 
alytically as factorized product states involving nuclear 
eigenstates of the (inhomogeneous) lowering operator. In 
analogy to the homogeneous case, such solutions cease to 
exist after the critical point at which we find diverging 
nuclear squeezing. These results are supported by nu- 
merical simulations that confirm the analytical consider- 
ations and provide further indications that other features 
of the phase diagram aside from the second order transi- 
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tion can be found in the inhomogeneous model. 

In order to take into account inhomogeneities in the hy- 
pcrfine coupling, we replace the homogeneous spin oper- 
ators of Eq. Q with inhomogeneous operators I a — > A a 
(a = x,y,z). We approximate the actual distribution 
of coupling strengths by n shells of spins with identical 
coupling 



(61) 



where A& represent homogeneous spin operators within 
the ith shell. Each homogeneous shell is assumed to be 
in a symmetric subspace Jj. 

In analogy to the homogeneous case we can con- 
struct approximate eigenstates of the lowering operator 
A~ | a) = a|a). To this end we perform a Holstein- 
Primakoff transformation on the homogeneous spin op- 
erators within each shell and displace the respective 
bosonic mode hi by /3j and expand the resulting opera tors 
in orders of 1/ 



A2 



"Ji- As we demonstrate in Appendix 
the choice of a particular displacement j3i uniquely de- 
fines the squeezing of the respective mode 6j if we de- 
mand that the corresponding state is an A~ eigenstate 
to second order in the expansion parameters, i.e., of or- 
der 0(^2 i l/ Ji). The corresponding eigenvalue is then 
given as a = X)f =i 9iVkjflj {ki — 2~ |/3i| 2 ). As discussed 



in Section iHIBl \ip) = |!) ® \a) is a steady state of the 
evolution to second order, if a = J^i diVkifti = — Jf2/^o- 
In contrast to the homogeneous case (n = 1) the latter 
condition does not determine the steady state uniquely. 
Several sets of displacements within the different shells 
can fulfill the steady state condition. However, all these 
microscopic realizations lead to the same macroscopic be- 
havior of the system such as the locking of the electron 
inversion (S z ) — 0. Furthermore, at the critical point, 
the solution is unique again (/% = 1 for all shells) and 
the considerations on entanglement of Appendix | A l| can 
be straightforwardly generalized to the inhomogeneous 
case with the result that also here at the critical point 
the entanglement in the system diverges indicating a sec- 
ond order phase transition. Obviously, above the critical 
point no such solution can be constructed and the system 
observables change non-analytically. 

Fig. [6] shows numerical results which confirm the above 
considerations. We find numerically the exact steady 
state solution for a model of two inhomogeneously cou- 
pled shells (gi = 2^) of size — 8 (broken lines), 
as well as for a system of 5 nuclear spins with coupling 
strengths ({sj l= i... 5 = {0.67,0.79,0.94,1.15,1.4}, dot- 
ted lines) . For low driving strengths Q we find the Over- 
hauser field building up linearly, as expected. The emer- 
gence of the thermodynamic phase transitions can be an- 
ticipated already for these low particle numbers. 

These analytical and numerical arguments for the 
emergence of a second order phase transition in the in- 
homogeneous case, suggest the possibility to find other 



features of the homogeneous phase diagram also in inho- 
mogeneous systems, such as NV centers in diamond. 

Another attractive realization of a central spin system 
is provided by singly charged semiconductor quantum 
dots: up to several 10 4 nuclear spins are coupled to a cen- 
tral spin-1/2 electron, driving and spin pumping of the 
electronic state have been demonstrated experimentally 
with high efficiency [29l|49]. In this setting, however, the 
inhomogeneity of the hyperfine coupling and the absence 
of an m = central spin state lead to a situation in which 
the effective nuclear Zeeman term Hi in Eq. becomes 
inhomogeneous (it is composed of Knight field, nuclear 
Zeeman energy, and the (homogeneous) detuning) and 
does not vanish for any choice of parameters. Therefore 
the above argument for a persistence of the second or- 
der phase transition does not apply. However, critical 
phenomena similar to the ones described above were ob- 
served in optically driven quantum dots |24) . The adap- 
tation of our model this and other more general settings 
is subject to future studies. 



VII. CONCLUSIONS 

In analogy to closed systems where critical phenom- 
ena arise from non-analyticities of the Hamiltonian low 
energy spectrum, in open systems critical phenomena 
are intimately related to the low excitation spectrum of 
the Liouville operator. We investigated a generic driven 
and damped central spin model and its rich steady state 
behavior, including critical effects such as bistabilities, 
first and second order phase transitions and altered spin 
pumping dynamics. We developed a two-step pertur- 
bative theory involving the expansion of nuclear fluctu- 
ations up to second order in a self-consistent Holstein- 
Primakoff transformation and the subsequent adiabatic 
elimination of the electron degrees of freedom in the 
vicinity to the steady state, which enabled us to provide 
a complete picture of the system's phase diagram. Link- 
ing common ideas from closed system phase transitions 
to the dissipative scenario, we were able to introduce a 
classification of the different transitions in the phase di- 
agram. 

The relevance of the considered model involves two as- 
pects. On the one hand Eq. ([I]) describes a simple yet rich 
model, which displays a large variety of critical phenom- 
ena. The limitation to symmetric states allows for an ef- 
ficient (and in the thermodynamic limit exact) perturba- 
tive treatment that gives deep insights into the nature of 
dissipative critical phenomena from a fundamental point 
of view. On the other hand the central spin model is 
general enough to have realizations in a large variety of 
physical systems (e.g., quantum dots, Nitrogen- Vacancy 
centers) . Our understanding of the critical phenomena in 
this model could provide insight into recent observation 
of critical behavior in such systems [24| [25]. Further- 
more the main features of the phase diagram discussed 
above can also be found if the central (two-level) spin is 
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replaced by a different physical system, e.g., a larger spin 
or a bosonic mode. The theory developed in Section IV 
can straightforwardly be adapted to different scenarios 
and opens the possibility to study dissipative critical ef- 
fects in a variety of different physical systems [13] ■ 

Finally, we showed that in a more realistic adaptation 
of the model incorporating an inhomogeneous hypcrfine 
coupling, the second order phase transition persists, in- 
dicating the possibility that the phase diagram remains 
qualitatively correct in this experimentally more realistic 
case. A more thorough analysis of the effects of inhomo- 
geneities is subject to future work. 



Acknowledgments 

We acknowledge support by the DFG within SFB 631 
and the Cluster of Excellence NIM (E. M. K., G. G. 
and I. C), the NSF (M. L. and S. F. Y.), CUA and the 
Packard Foundation (M. L.), as well as the ECR (A. I.) 
and AFOSR under the MURI award FA9550-09-1-0588 
(S. F. Y.). 



Appendix A: Approximate Eigenstates of the 
Lowering Operator 



1. Homogeneous case 



In Section IIII Bl we have seen that we can construct 
the exact steady state along segment x if we assume the 
nuclear system to be in an eigenstate of the spin lower- 
ing operator I~ \a) — a\a). Although it readily can be 
show the this operator exactly features only the eigen- 
value a = 0, we can construct approximate eigenvalues 
in an expansion in 1/ J. 

We stress the point that in the bosonic analogue eigen- 
states of the annihilation operator are coherent mini- 
mum uncertainty states that display no squeezing. As 
we will see, the eigenvectors of the atomic lowering oper- 
ator in contrast are squeezed coherent atomic states (on 
the southern hemisphere of the Bloch sphere) , where the 
squeezing parameter depends uniquely on the rotation 
angle of the Bloch vector. 

As noted in Section IIVI the Holstein-Primakoff trans- 
formation Eq. Q provides an exact mapping between 
spin operators and a bosonic operator in the subspace of 
total spin quantum number J. In the following we show 
that approximate eigenstates of the lowering operator I~ 
can be expressed as a squeezed and displaced vacuum of 
the bosonic mode b 



where D(f3) 



D(fi)S{-r(J3)) |0) =: , (Al) 



are the displacement and squeezing operators, respec- 
tively and |0) = \J — J) the fully polarized nuclear state. 



We find the squeezing parameter uniquely defined by the 
displacement r = r(/3). 

Without loss of generality we assume /? £ R (and thus 
r £ R), i.e., the Bloch vector lies in the x — z plane. 
General states (3 £ C with arbitrary Bloch vectors on the 
southern hemisphere, can straightforwardly be derived by 
a rotation around the z-axis. Note that the correspond- 
ing states on the northern hemisphere can be constructed 
accordingly as eigenstates of the ascending operator I + . 



In order to show that Eq. ( Al ) defines an approximate 
eigenstate of I~ we first consider the transformation of 
the nuclear operator under the displacement and squeez- 
ing operator. Recall that according to Eq. ([9| the dis- 
placed nuclear operators can be expanded in orders of 

D\p)I~D{p) (A2) 
2 J - (&t + y/j/3*)(b + VI j3) (b + \fjp\ 



where 



= J Jo +VJj^ + G{l), 



J ~ = VkP, 

Jf = y/2(l -P*)(jib + Ub*) 

= y/2(l- /3*)&(r)bS(r), 

?*—S? — an d sinh(r) 



(A3) 
(A4) 



and cosh(r) 

/3 2 

2 x /2fc(l-/32)' 

to complex /3 is straightforward and leads to Eq. (|44j)). 
Thus it follows 



2V , 2fc(l-/3 2 ) 

which defines r = r{0) (the generalization 



S*(-r)Di(fi)rD(P)S(-r)\0) 

= JJ -\0)+O(l), 



(A5) 



since b |0) = 0. 

Multiplying both sides by D(/3)S(—r) yields the de- 
sired approximate eigenvalue equation 



r 



0(1) 



(A6) 



In the thermodynamic limit the term 0(1) is negligible 
and the eigenvalue equation is exact |57j . 

Using the above representation we study the spin prop- 
erties of the states \a). In the following all expectation 
values are understood to be evaluated in the squeezed 
coherent state \j3): (O) = {0\O\P). 

Straightforwardly, one derives the nuclear mean polar- 
izations 



(4) 

a,) 



i 

2 
1 

2i 

J{P 2 



(3\(I + +I-)\(3) 

w\(i + -n\/3 
) + 0(l), 



jVkp + 0(l), 
= + 0(1), 



(A7) 

(A8) 
(A9) 
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where in the last equation we used the expansion 



Eq. (14). Note that the Bloch vector is orthogonal (up 



to order 0(1)) to the y-direction for all (real) a and of 

length | (I) | = \/(^) 2 + (I y ) 2 + (h) 2 = J + 0(1)- 

Using Eq. ( A6 1 and the angular momentum commuta- 



tion relations one readily calculates 



1 



M) = —Ah) +0(1), 



infinitely many microscopic realizations (i.e. sets of Bi) 
that fulfill a = JiVkiPi- Only the maximal eigen- 
value a = J features a unique steady state that displays 
diverging squeezing as one readily shows analogous to the 
homogeneous case. 



(A10) Appendix B: Rotated Squeezed Thermal Spin States 



= 2./(i 0(1), 



Wl-(V^9) 2 + 0(l) 



A key concept of the paper are RSTSS, a generaliza- 
tion of squeezed coherent spin states to mixed states, 
parametrized via an effective temperature. They describe 
nuclear states which are fully polarized, rotated, and fea- 
ture fluctuations which can be described by a bosonic 
mode in a thermal (potentially squeezed) Gaussian state. 
In Section |IV A| we show that the truncation of every 
nuclear operator to a subspace of total spin J can be 
expressed in terms of a bosonic mode b and its displace- 
/ ~ ment B € C, using a Holstein-Primakoff transformation 

e y = 2(AI 2 )/\(I)\ = yj\ {VkPf + 0(1/ J). (All) [compare Eq . j| j Eq . Q] 



where as usual (AO 2 ) :— (O 2 ) — (O) 2 and we used the 
identity 1 - (Vk(3) 2 = (1 - B 2 ) 2 . 

Thus, we find for the squeezing parameter in y- 
direction, 



The squeezing diverges for the state that realizes the 
maximal eigenvalue of the lowering operator {ykB = 1). 
This corresponds to a state fully polarized in x direction. 



2. Inhomogeneous case 

We approximate a system of inhomogeneous hyper- 
fine coupling by grouping the nuclear spins into n shells. 
Within a shell I the nuclear spins have identical coupling 
cji and the respective (homogeneous) spin operators Aa 
(a = x,y,z) are truncated to a symmetric subspace Jj. 
The total spin operators can then be written as 



(A12) 



We define collective displacement and squeezing opera- 
tors 



V = n 



n ps/Jifiibl-s/Ji^bi 



(A13) 
(A14) 



where the bi is the respective bosonic operator for shell i. 
Also here the squeezing parameter r» depends uniquely 
(with the same functional dependence as before, cf. 
Eq. |A4| )) on the displacement Bi within the shell, if we 
demand the first order in the eigenvalue equation to van- 
ish: 

A-VS |0> = {J2 JiVkiPi)VS |0) + 0(1), (A15) 



where ki = ^2 — B 2 and |0) = |0)® n is the vacuum of 
the shell modes. 

We emphasize that in general the eigenvalues are 
highly degenerate. For a given eigenvalue a there are 



(Bl) 



where e = l/\/J, and the bosonic operators J™ con- 
tain combinations of products of n bosonic operators 
b, b* . Jq € C, describes the semiclassical expectation 
value which is fully determined by the displacement B. B 
quantifies a rotation of the fully polarized nuclear state 
on the Bloch sphere. The higher order operators J£ 
(n > 0) describe quantum fluctuations around this semi- 
classical nuclear state. RSTSS are those states where the 
mode b is in an undisplaced ((b) = 0), squeezed thermal 
state, which is fully determined by its covariance matrix 
r [Eq. (41)]. These bosonic states constitute the natural 



steady states of the quadratic Master Eq. (25), and we 



find in Section IV B that across the whole phase diagram 
one steady state of the system can always be described 
as a RSTSS. 

Note that in the limit where the effective tempera- 
ture of the Gaussian state is zero, we recover the class of 
squeezed coherent spin states [33], which constitute the 
solution along segment x. 



Appendix C: Solving Eq. (24| 



In order to find the solutions to Eq. (24) (which are 
numerically difficult to find) we first note that 

(A) ss = 0^(b) = (6t) = ^ (J,-) = = 0, (CI) 

where the time derivative is understood with respect to 
the first order Liouvillian 



C XP = - i[a(S x J? + S y J?) + (aS+S- + 8u)Jf,p], 



(C2) 
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and in the usual way we define 



1 



(C3) 
(C4) 



Using the relation J{,J{ — if-ijkJo one finds the equa- 
tions 

= (Jf ) = a ((S y ) ss J Q z - (S z ) ss J y ) - u:jy, (C5) 
= (jf) = -a ({S a )„J$ - (S z )ssJ Q x ) + ujJ x , (C6) 
= (Jf > = a {(Sy)„J$ - (S X ) SS J$) . (C7) 

Furthermore from the definitions of the 's one finds 

l = (J?) 2 + (J?) 2 + (Jo z ) 2 - (C8) 

The steady state expectation values (Si) ss are found 
directly via [cf. Eq. Q] 



CaP=l{S pS + ~-{S + S ,p}+) 



(C9) 



- *[S x (2n + oj?) + o5 v .7 » + aS+S-Jf, p], 
by solving the resulting optical Bloch equations 



(CIO) 



= -^ x )+a^(5 z )-aJ ^), 

= -|(5„> - (2n + aJ x )<5 2 } + aJo z <S !C ), (Cll) 

= -7((S,) + 1/2) + (2fi + aj 2: )(^) - aJ^S*). 

(C12) 

This set of coupled Bloch equations for the six vari- 
ables {(Si), Jl } can be solved analytically. The solutions 
which feature second order stability (see Section IV B 1 ) 
are displayed in Fig. [2] Via Eq. ( |l"o| and Eq. ( 14 1 can 
be deduced unambiguously from a given set {{Si}, J I }. 



Appendix D: Deriving t he s econd order term of 
Eq. ([20 » 



The first term of the second order of Eq. ( 20 ) is of the 



same form as the first order and can readily be calculated: 

Tr s (PC 2 Pp) = -i[a/2((S+) ss Ji + (S-) SS J+) (Dl) 
+ (a(S+S-) ss + Suj)J 2 z ,a], 
= -i[B*b 2 + B{tf) 2 + Ftfb, a], 

with the /3-dependent coefficients (remember that also 
the electron steady state expectation values are functions 
of 0) 



B = — 
F = — 



(10 



16Vfc 3 
a 



[(4k + \0\ 2 )(S-)ss+P 2 (S+), 



(4k+\l3\ 2 )(0(S+) ss +0*(S-) t 



(D2) 
(D3) 



+ a({S+S-) ss + 6u/a). 



Next, we consider the second term of the second order 
perturbative master equation 

-Tr^PdQC^QdPp) (D4) 
= - Tr s (P£i(l - P)Lq \l - P)C x Pp) 

= / dTTr s (P£ ie CoT XiPp) 
Jo 

/>oo 

- / dTTr^PdPCtPp), 
Jo 

where we used the Laplace transform — Cq 1 = J °° dTe c ° T 
and the property e c ° T P — Pe c ° T = P. 
Noting that 

Tr s (PdX) = -iTr s ([bA + tfA\X}), (D5) 
and using Eq. ( p2| we find 

/>oo 

- / dTTr.(PCiPCiPp) (D6) 
Jo 

dT(A a ) ss (A?) ss [b a ,[bP,a\] , 



where a, = {, 'void' and the Einstein sum convention is 
used. 

In the same fashion we find 

/>oo 

/ dTTr a {PC x e c ° r L x Pp) (D7) 
Jo 

dr(A a (r)A^(0)) ss [b a , [bP , a]] 

dr([A a (T),A^(0)]) ss [b a ,ab?]. 

Here we defined the autocorrelation func- 
tions (A°'(t)A' 3 (0)} ss = Tr s (A a e CoT Af } p ss ) and 
([A a {T),Af>(0j\) aa = Tr s (A a e c ^[A^ Pss ]) (cf. e.g., 
\5U\ pp. 22). 

Putting together the results Eq. (|D4| reduces to 



-Tr^PdQC^QdPp) 



(D8) 



dr(AA a (T)AA^(0)) ss [b a , [6* a]] 
dr([AA a (r),AA^0)]) ss [b a ,ab ] 



AO := O — (0) ss . Since we choose the displacement 
such that (A a ) ss = [Eq. ([22b] it is AA a = A a . Merging 



Eq. (|D1[) and Eq. (D8), and regrouping the terms, one 



readily derives equation Eq. ( 25 1 



Calculation of the coefficients 



In order to determine the coefficients Eq. ( [28] ) we have 
to calculate terms of the kind / °° dr{AA a (TjAA p (0)) aa 



22 



and / o °°dr(AA a (0)A^(r)) ss . Exemplarily we will cal- 
culate the two terms for a — (3 — 'void'. 

First, defining v = (^=(2fc - \P\ 2 ),-^,(3a) T 

we can write AA = v* ■ AS (and with w = 
(-iTfc^*) 2 ' ^(2fc-|/3| 2 ),/3*a) T wefind AA* = w*-AS 
). Likewise it is AA^ = A& ■ v (A A = AS^ ■ w). 
Consequently we compute: 

/>oo 

/ dr(AA T AA) ss (D9) 
Jo 

=v* (J dT(AS T AS^) ss ^j w 

=v* uf°° dTe MT (ASA&) 8S ^ w 
=v* (-M^iASAS^)^ w = tT^itZT, 
where we applied the Quantum Regression Theorem in 



we write 



the second step and used the definitions of Section ( IV A I . 
Noting that 

roo _ / roo \ t 

j dr(ASASi) ss = I J dT(AS T A^) ss ) (DIO) 



- (ASA&) SS M^ =T 2 = Tl 



roc 

/ dr(AAAA T ) ss = v*T 2 w. (Dll) 
Jo 



Analogously, we find the relations 



/•DO 

/ dr(AAtAA) ss -- 
Jo 


= w* Fiw 


/>oo 

/ dr(AA^AA T ) ss -- 
Jo 


= W* F2W 


POO 

/ dT(AA T AA^) ss -- 
'o 


= v* Tiv, 


poo 

/ dr(AAAA\) ss 
Jo 


= V f 2 V, 



(DI2) 



such that all coefficients of the effective Master Eq. ( 20 ) 
can be calculated by simple matrix multiplication. 
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